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PREFACE 


The report contains the results of several papers 
related to modeling using a class of the bivariate gamma 
distribution. The separate papers contain loosely related 
subjects pertaining to this problem. Since the separate 
papers were prepared at different times during the contract 
period and have been submitted for publication in the open 
literature and each paper is intended to be self-contained, 
there is some redundancy in tables and illustrations. 

Each of the papers in this report were extensions 
and/or generalizations of the results given in NASA TM- 
82483, entitled "A Bivariate Gamma Probability Distribution 
with Application to Gust Modeling," by 0. E. Smith, S. I. 
Adelfang, and J. D. Tubbs. A modification of this paper is 
currently under review by Communications in Statistics . 

The first paper in this report, entitled "A Note on 
the Ratio of Positively Correlated Gamma Variates," has 
been accepted for publication in Communications in Statis - 
tics and it presents some new analytical results using a 
class of the bivariate gamma distribution. Comparable 
results were available in the open literature using a 
different cxass of the bivariate gamma. 

The second paper is entitled "A Method for Determin- 
ing if Unequal Shape Parameters are Necessary in a Bivariate 
Gamma Distribution" and is an application of the results 
given in the first paper and addresses questions concerning 

iii 



hypothesis tests for equality of shape parameters from 
correlated gamma distributed variates. This paper is 
currently under review by Teclinometrics . 

The third paper, entitled "A Differential Equations 
Approach to the Modal Location for a Family of Bivariate 
Gamma Distribution," contains extensive analytical results 
for the location of the mode as a function of the free para- 
meters. To the authors' knowledge this is the only such 
representation for a non-gaussian bivariate distribution. 
This paper has been submitted to SIAM J . on Scientific 
and Statistical Computing . 

The fourth paper is a report summarizing the analysis 
of some wind gust data using the analytical results devel- 
oped in relationship to the modeling application. 
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CHAPTER I 


A NOTE ON THE RATIO OF POSITIVELY CORRELATED GAMMA VARIATES 


J. D. Tubbs 

Department of Mathematical Sciences 
University of Arkansas 
Fayetteville, Arkansas 


0. E. Smith 


Systems Dynamics Laboratory 
NASA Marshall Space Flight Center 
Huntsville, Alabama 


ABSTRACT 

Mielke and Flueck (1976) derived the density function and 
corresponding moments for the ratio of correlated gamma distrib- 
uted variates. They considered a class of bivariate gamma dis- 
tributions suggested by Cherian (1941) and David and Fix (1961). 
Recently, Lee, Holland, and Flueck (1979) derived some additional 
distributional results using this class of functions. This paper 
derives similar results using a different class of bivariate 
gamma distributions. 


1. INTRODUCTION 

Mielke and Flueck (1976) derived the distributional results 
for the ratio, R, of correlated gamma distributed variables. 

There are several classes of the bivariate gamma distribution [three 
are summarized in Mardla (1970) and an additional two in Johnson and 



Kotz (1972)]. Mielke and Flueck (1976) derived the distributional 
results for the ratio, R, of correlated gamma distributed variables 
using the Cherian-David-Flx class of bivariate gamma random variables 
[Cherian (1941) and David and Fix (1961)]. That is, let X, Y, and P 
denote Independent gamma random variables with common scale parameter 
X and respective shape parameters a - 3 “ and for 0 < ^ < 

min (a ,g) . 

Then it can be shown that the bivariate probability density 
function for U X + P and V = Y + P is given by 

min (u,v) 

-_r_1 o_r_l 

dp (1.1) 


J - min (u,v) 

5-1, .a- 5 - 1 , V0-5-1 p 

0 p (^"P) ® 


for K = r(o-5)r(e-5)r(5) 


when the scale parameter X is assumed to be unity. Mielke and 
Flueck (1976) showed that (1.1) can be written as 




a-1 e-5-1 -(u+v) 
i V e 

r(a)r(B- 5 ) 


a- 5-1 B-1 -(u+v) 
u V e 

r(o-5)r(B) 


F*(5 , 1 + 5 -a :B.u/v,-u) if 0<.u<v 


F*(5 , 14^-a : B»v/u,-v) if 0<v<u 


( 1 . 2 ) 


OB m n 

» / . \ _ n I m n id n , n ^ 

where F, (a,b,c:x,y) » E 7 —;: ; — r x y, is a degenerate 

m,n =0 m+n 

two variable hypergeometric function [Gradshteyn and Ryzhik 
(1967), p. 1067] and “ r(a+n)/r(a). Thus, U and V are gamma 

random variables with shape parameters a and B and positive depend- 
ence parameter 5. In particular, E(U) = Var(U) = o, E(V) * Var(V) 

” B, and Cov(U,V) = 5 . 

Mielke and Flueck (1976) derived the density function for 
R =■ U/V using a change of variables. That is. 
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fR(r) 


r°~^(l+r)^~° ^ F (C,a+8-C,l4^-a,a:r/l+r,r) 
B(a-C.a) 

if 0 < r < 1 


(1.3) 


B(a-£,8) 


Fj(C,a+6-5, l+5-a,3:l/l+r,l/r) 
if 1 < r 


* (a)_j_(b) (c) 

where F^ (a,b,c,d;x,y) = Z — — x y , |x| < 1 |y| < 1 

m,n*0 ^ ^m+n 

is a two variable hypergeometric function [Gradshteyn and Ryzhik, 
(1967), p. 1053], and B(a,b) = F (a)T (b)/F (a+b) . In addition, they 
show that the integral moments of R are given by 


E(R^) - jl, (p 




(e-j) 


for s > 0. 


(1.4) 


s 


In particular, 

E(R) = ag-5/B(B-l), B > 1 


(1.5) 


E(R^) 


g(C+l) . 2(o-5)£ (a-C+l)(g-£) 

8 ( 8 + 1 ) 8 ( 6 - 1 ) ( 8 - 1 ) ( 8 - 2 ) ’ 


I’^icently, Lee, Holland, and Flueck (1979) were able to obtain 
comparable results for density of R using the Cherian-David-Fix 
class of densities by expressing f_ as a weighted difference of 
hypergeometric functions. The purpose of this paper is to derive 
comparable results for R using a different class of the bivariate 
gamma distribution. This class is a special case of the one 
suggested by Jensen (1970) as modified bv Gunst and Web;.;ter (1973), 
The next section contains a brief discussion of this class of 
distributions. In section 3 the derivation of f is given using 
this class of functions. Section 4 outlines a possible applica- 
tion for the probability function in the area of hypothesis 
testing for the equality of shape parameters in the presence of 
correlation. 
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2. 


GUN5T AND WEBER CLASS OF BIVARIATE GAMMAS 


Gunst and Weber (1973) proposed a computationally feasible 
method for deriving the joint density function for the bivariate 
chi-square distribution. Since the chi-square is a special case 
of the gamma, this method was used for the bivariate gamma case. 
That is, a bivariate gamma density function for U and V with 
common scale parameter X = 1 and shape parameters a, 3, (a < 3) 
is given by 


f(u,v) 


^a-1^3-lg-[ (u+v)/(l-n)] “ “ 

(i-n)“r(a)r(3^) j k 


(3-g+k) (uv)'^v^ 
(l-n)^^'*'^(3+j+k)j!k! 


( 2 . 1 ) 


where n = p ^ (3/a) , P is the correlation coefficient between the 
variables U and V. Gunst and Webster (1973) suggested this class 
of densities in that they are computationally tractable and do 
net involve mathematical functions, such as Laguerre polynomials 
or convoluted sums [Jensen (1970) and Kibble (1941)]. Smith, 
Adelfang, and Tubbs (1982) discuss this class of densities in 
greater detail. 

In the next section the distributional properties for the 
ratio, R, are derived using the Gunst-Webst ?r class of bivariate 
gammas . 

3 . RATIO OF CORRELATED GAMMA VARIATES 

By letting R = U/V and S = U+V, the joint pdf for R and S 
can easily be shown to be 




QO OD 




E I c. 

j k 


(i«)^ 



a+j-1 



3+j+k-l 

e 


-s/(l-ri) 


(3.1) 


where Cj^ = [ (l-n)°'r (o) (3-<i) ] c„ = 


n^'*'^(3-a-t-k) 


^ (l-n)^^'^'V(3+j+k)j!k! 

by Integrating over S the pdf for R becomes 


Hence, 
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00 


(3.2) 


( 1 )* 


f (r) - (1-n)® I E c c 

^ a_A 1 n j ^ 


j“0 k“0 


where =j - (a)jn^/jl, Cj^ » (B-a)j^nVk!, (a)^ - r(a+n)/r(a), and 
B(a,b) - r(a)r(b)/r(a+b). 

Whenever the shape parameters are equal then the de.isify func- 
tion for R is given by 


fj^(r) = (1-n)® E c B(a+j 


(3.3) 


j=0 


From (3.2) it can be shown that the m^^ raw moment for R is 
given by 

00 00 

E(R™) =* (1-n)^ E c E c. B(a+j+m,S+j+k-m)/B(a+j ,6+j+k) (3.4) 
j-0 ^ k-0 ^ 

if m < B. In which case, it follows that 


00 GO 


E(R) - (1-n)^ E c, E c. (ci+j)/(B+j+k-l) 
j ^ k 


(3.5) 


and 


OO 00 


E(R^) - (1-n)^ E c I c^ (a+j)(a+j+l)/(B+j+k-l)(B+j+k-2) (3.(. 

j ^ k 


Whenever n “ 0, then 


E(R) = a/(6-l), E(R ) = a(a+l)/(B-l)(3-2) 


(3.7) 


which agrees with the values given by Mielke and Flueck (1976) 
whenever 5 “ 0 ana with Lee, Holland, and Flueck (1979) whenever 
a - 0. 

Lee, Holland, and Flueck (1979) discuss some of the mathe- 
matical properties for the density of R for various values of 
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a,0, and n. They demonstrated that the density can be «> at r"l 
whenever either of the shape parameters is less than one. How- 
ever, in the Gunst-Webster construction by assuminf that a > 1 
and a < g the density fvinction given in equation (3.2) is stable. 
Figures 1-4 Illustrate the various shapes that f„(r) has as a 
function of the three parameters. 

A definite computational advantage of equation (3.2) versus 
equation (1.3) stems from the ability to compute the tail proba- 
bilities for R. By letting a=a+j and b=g+j+k, we have 

OD 00 

Fg(ro> - a-n)» £ Plr2.,2b S- < 3 .s) 

j k 

where F denotes a random variable from an F-distribution with 
r,s 

r and s degrees of freedom. Note if n = 0» then (3.7) becomes 


which agrees with the well known results concerning the ratio of 
independent chi-squares. Furthermore, if n 0 and a = g then 
(3.7) becomes 




(1-n) r c P[F„, . , - < r^] 

j j 2va+j),2(a+j) = O' 


(3,10) 


which is similar to an expression given by Johnson and Kotz (1972'), 
Chapter 40, Section 3. 

4. APPLICATION 

In this section an application is given for computing the 
cdf of R, given by equation (3.7). Diagram 1 defines the area 
given in equation (4.1). 
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r • u/v 



By letting cot 8 q ■ U/V » and G(9) - l-Fj^(rQ), one has 


G(8) - (1-n) 


8 


j ^^^2b,2a “ 


(a/b)tan 8] 


(4.1) 


Figure 5 contains the graph of the function G(8) versus 8 for a ■ 1 
and 8 ■ 1, 2, or 3 and n “ 0, .25, .50, and .75. From this figure 
and other cases which are not Included one observes that whenever 
0=8 then G(A5°) = .5 and G(45°) < .5 whenever o < 6. This obser- 
vation and addi clonal properties were used In developing a test for 
the hypothesis 


Hq: o = 8 vs. o < 8 


(4.2) 


The procedure Is presented In Tubbs (1983) and uses the Cramer- 
Von Mlses criteria for testing (4.2). That is, define 

= n y* {F^(r) - F^(r)}^dFj^(r) (4.3) 

where Fg^(r) is the cdf for the null distribution given in (3.10). 

F (r) is the empirical distribution for r. = u. /v. and the r. *s 
are arranged in Increasing order. Whenever Is true, then is 
distribution free and has a convenient computational form given by 


W 

n 


1 

12n 


n 

+ Z {Z. 
i=l 


(2i-l)j ^ 
2n ^ 


(4.4) 


where Z. = F„(r. ). is rejected if W exceeds a specified 
X H X u n 

critical point. Tubbs (1983) considers the properties of this 
test procedure in greater detail. 


5. CONCLUSIONS AND SUMMARY 

This paper derives both the density and the distribution func- 
tions for the ratio of positively correlated gamma variates using 
a modification of Jensen's bivariate gamma distribution. The 
expression for the moments differ from those given by either 
Mielke and Flueck (1976) or Lee, Holland, and Flueck (1979). 
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However, all the expressions are identical whenever the variates are 
uncorrelated. A principal advantage found in this representation 
stems from the ability to compute the CDF of the ratio. The value of 
the CDF for the ratio was shown to have potential application to the 
problem of testing for equality of shape parameters in a particular 
family of the bivariate gamma distribution. 



Diagram 1 
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CHAPTER II 


A METHOD FOR DETERMINING IF UNEQUAL SHAPE 
PARAMETERS ARE NECESSARY IN A BIVARIATE 
GAMMA DISTRIBUTION 


J. D. TUBBS 

Department of Mathematics 
University of Arkansas 
Fayetteville, Arkansas 


ABSl "ACT 

A procedure for aiding an experimentalist in decid- 
ing between four and five parameters in a Jensen's type 
bivariate gamma distribution is presented. The procedure 
is based upon the properties of the CDF for the ratio of 
correlated gamma distributed variates. The criteria of 
interest is posed in a test of hypothesis setting and 
results are presented using the Cramer-Von Mises test 
of fit. 

1. INTRODUCTION 

Smith and Adelfang (1981) discuss the applicability 
of a bivariate gamma distribution as a parametric model 
for wind gust amplitude and length. In modeling this 
bivariate data with a gamma distribution, it was neces- 
sary to find a distribution that would allow for correla- 
tion between the random variables X and Y when the marginal 
distributions are univariate gammas with possibly unequal 
shape and scale parameters. That is, X ~ and 



Y “ G(y,,, 3 ) where the probability density function for 
j y 

Z ~ G (y , 8 ) is given by 

f^Cz) = 6^2^"^e"^^/r(Y). (1.1) 

A brief survey of the open literature reveals that 
there are several classes of the bivariate gamma distribu- 
tion. One need only consult Mardia (1970) and Johnson and 
Kotz (1972) to find five classes of the bivariate gamma 
distribution [Kibble (1941), Cherion (1941), McKay (1934), 
Jensen (1970), and Moran (1969)]. Of these classes only 
Jensen (1970) and Moran (1969) allow for unequal shape 
parameters and both of these have computational limitations 
which affect their utility to the experimentalist. Recently, 
McAllister, Lee, and Holland (1981) and McAllister (1983) 
have addressed the limitations with Jensen's model and 
provided results which overcome many of the computational 
difficulties. However, at the time of Smith et (1983) 
development these results were not available. Hence, they 
modified a bivariate chi-square model given by Gunst and 
Webster (1973), This allows for possibly unequal shape 
parameters and is computationally tractable. The model is 
not as general as that given by Jensen (1970) , however, one 
can derive the bivariate model given by Kibble (1941) as 
a special case whenever the shape parameters are equal. In 
this paper, the unequal shape parameter model will be 
referred to as the five-parameter model and the equal shape 
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case as the four-parameter model. Smith, Adelfang, and 
Tubbs (1983) discuss the properties of these distributions 
and it is apparent that the four-parameter has numerous 
computational advantages over the five-parameter model. 

So if one assumes that the data is correctly modeled by 
this class of the bivariate gamma distribution, a question 
of practical interest becomes. How does one decide if the 
five-parameters are really necessary? The purpose of this 
paper is to present a procedure which would aid the exper- 
imentalist in answering the above question. The problem 
is posed in a hypothesis testing setting. That is, test 
the hypothesis 



versus 



^ ^y 


(1.3) 


It should be noted that the proposed method is not an 
omnibus test of fit for the bivariate gamma against all 
other possible models. Instead the procedure is intended 
for deciding between the four or five parameter models as 
given in Smith, Adelfang, and Tubbs (1983) . 

The next section contaixis the distributional results 
needed for the test of hypothesis (1.2). The test proce- 
dure is given in section 3 and evaluated in section 4. 
Section 5 contains a summary and remarks concerning some 
of the limitations of the procedure. 
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2 . 


DISTRIBUTIONAL RESULTS 



Smith, Adelfang, and Tubbs (1983) modified a bivariate- 
Chi square distribution given by Gunst and Webster (1973) 
and obtained the density function given by 

f(x,y) = K,/K^ t f c., x'^(ny)^'‘'^ (2.1) 

^ ^ j=0 k=0 

where 

Yv"! Y„-l 

* X * y y exp-{ (x+y)/ (l-n)> , 

K 2 = (1-n) (Yj^)r (Yy-Y^) * 


and X = XB , y ■ YB^, 8„, are known scale parameters, 
X y X y 


n = p / Y /y » P is the correlation coefficient between the 
y X 

variables X and Y. The joint probability distribution func- 
tion is given by 


F(x^,y^) = Pr(X < x„, Y < yj 

= J I E d,, H(y +j,x /(1-ri)) 

j=0 k=0 

•H(Yy+j+k, y^/d-n)) (2.2) 

where 

J = (1-n) ^/r (Y^)r (Yy-Yjj) . 

d-y. = n^'^^r(Yy-Yx+k)/r(Yy+j+k) jl kl 
H(a,x) = e“^ dt. 
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Equations (2.1) and (2.2) are for the unequal shape para- 
meters and will be referred to as the five-parameter model. 

It should be re-emphasized that this model is not completely 
general in that one assumes that y > y„ and the correlation 
between variables X and Y are restricted to the interval 


[0, ti/yITy^] for ne[0,l]. 

A y 


If Y = Y 
'x 'y 


Y then it can be shown that (2.1) and 
(2.2) reduce to the well known functions g,iven by Kibble 
(1941). That is, the density function is given by 
f(x,y) = (xy)^“^exp-{ (x-^y)/(l-n)}/r (y) 


• f (nxy)/(l-n)^)-^/r(Y+j) jl 
j=0 

and the distribution fmction becomes 

F(x,y) * (l-n)^/r(Y) r n^/r(Y+j) jl 

j=0 

•H(y+J |X/ (1-n)) H(y+J ,y/ (1-n) ) . 


(2.3) 


(2.4) 


Equations (2.3) and (2.4) will be referred to as 
the four parameter model. A comparison of the distribution 
function given in (2.2) and (2.4) reveals that there are 
distinct differences in terms of the computational com- 
plexity. Thus for computational reasons the experimental- 
ist would like to know how much greater does Yy have to 
exceed y^ before equation (2.2) is really necessary. 

Ideally he would like to answer this question before using 
both (2.2) and (2.4) then selecting the results which are 
more gratifying. In order to eJdress this issue, this 


II-5 



paper considers the problem of testing hypothesis (1.2) 
versus (1.3) using an univariate random variate given 
by the ratio of X to Y, R = X/Y. Tubbs and Smith (1983) 
derive the density and distribution functions for R when- 
ever the bivariate density is either (2.1) or (2.3). That 
is, if equation (2.1) holds then the density function for 
R is given by 

f„(r) » (1-T,)^y E E c.c,B(a,b)"V"V(l+r)^'^^ 

^ j=0 k=0 J 

where B(a,b) = r (a)r (b)/r (a+b) , Cj - (a)jn^/jl, 

Cj^ = (b-a)^n^/lcl, a = Yx+J» ^ = Yy+j+k, and (a)^ = r(a+n)/r(a). 
The distribution function for R is given by 


Fj^(r^) = Pr [X/Y < r^ J 


00 CO 


where F denotes a random variable from an F- distribution 
r ^ s 

with r and s degrees of freedom. The corresponding functions 


whenever Yv “ Y« “ Y given by 
X y 

fp(r) = (1-n)^ r c. B(a,a)“ V (l+r)^® 

^ j=0 J 

and 


(2.7) 


° jig ‘j f’^l''2a.2a i ’"q) 


( 2 . 8 ) 


where a = y + j • 



3. 


HYPOTHESIS TESTING 




Since R = X/Y is a mivariate random variable it is 

informative to graph versus r. However, since r > 0 

a more meaningful graph can be produced by letting 6 = 

cot ^r and G(e ) = 1-F, ) where 0^ “ cot . The area 

o' 1 o' o o 

corresponding to is shown in diagram 1. Further- 

more, it follows that 

Y Y “ OD 

G(e ) «= (1-n) ■'I Z c.c, Pr[F« ^ (a/b) tane J 
j=0 k=0 ^ ^ /a, ZD o 

(3.1) 

in the five-parameter model and 


G(e^) - 2^ < tanej 


(3.2) 


in the four-parameter case. 

Since 0 is restricted to the finite interval (0, tt/ 2), 
it is somewhat instructive to plot G( 0 ) versus e as func- 
tions of the free parameters, y , y and As in Tubbs 

X y 

and Smith (1983) the scale parameters are assumed to be 
known and hence equal to one. This restriction will be 
addressed later in the paper. Figures 1-3 contain some of 
the illustrative cases. 

From these plots one observes that G(45®) * .5 when- 
ever the four-parameter model holds and G(45°) < .5 in 
the five-parameter models. Rather than just using this 
observation a function was selected to measure the distance 
between these distribution functions. The Cramer-Von Mises 
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0 




o 















type goodness-of-fit procedure was selected since the test 
is distribution free whenever the parameters are specified. 
Furthermore the test statistic is easy to compute. 

Let 


TT / ^ 

'n “ - G^(e)}^dG(e) 


(3.3) 


where G(e) is given in (3.2), Gj^(9) is the empirical dis- 
tribution function of 0. = tan'^fr.), r.- = X. /Y. are 

r ■ X ^ ^ 

arranged in increasing order. Whenever hypothesis (1.2) is 
true, then has the convenient computational form given 
by 


n 


W = l/12n + Z {z. - 


r2i-l)^2 


n 


i=l 


TiT 


(3.4) 


where = G(e^). Furthermore, from Anderson and Darling 
(1951) one can reject (1.2) whenever W^ exceeds a specified 
critical point. These critical points are given from Ander- 
son and Darling's asymptotic distribution. Stephen (1976) 
defines a procedure for modifying the critical points for 
small samples, however, the underlying problem of modeling 
bivariate data will probably dictate large sample sizes. 


4. EVALUATION OF THE TEST PROCEDURE 


In this section the procedure defined in the previous 
section is evaluated. The evaluation is performed in two 
parts. The intent of the first part was to determine 
whether or not the procedure even works . That is , are the 
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apparent visual differences between the function G(0' as 
seen in Figures 1-3 significant in the "Cramer-Von Mises" 
metric. The second part of the evaluation concerns the 
robustness of the procedure to the nuisance parameters. 

In the first part, let 

7t/2 « 

D^(6) = n/ {G(0) - A( 0 )}^dG( 0 ) (4.1) 

o 

where G(0) is given in (3.2) and A(o) is given by (3.1) 
when Yv “ Y + <5 , for 5 >0. For positive integers n, 

y X 

compute . 

> D^(6)]. (4.2) 

If the alternative hypothesis given by 
“l- '<K ^ 'ly ~ ^ 

holds, then the expected value of In (3.3) is given by 
D^(d). Hence, a^(6) is the expected type I error of test- 
ing hypothesis (1.2) as a function of 6. Table 1 contains 
the value of “j^(<5) for various values of the parameters. 
The a^(6)'s were computed using Tiku's approximation to 
the asymptotic distribution of [Tiku (1965)]. 

For example, from Table 1 one would expect the test 
to reject integer (6=1) differences between the shapes for 
X and Y at the 95% significant level whenever n > 50. 

The procedure used to generate the values in Table 1 
is somewhat unconventional; however, they do indicate that 
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the test procedure would be sensitive to differences in 
the shape parameters that exceed unity. A Monte Carlo 
simulation was also performed. The results are not 
reported in the interest of space and since the simulation 
was quite limited. A detail simulation is very expensive 
due to the computational cust.^ in computing the null dis- 
tribution G(0) needed in evaluating type I errors. It 
xs especially costly to simulate any type II errors. In 
spite of these restrictions upon the simulation's merit, 
the results were supportive of the expected results given 
in Table 1. 

The second part of the evaluation is concerned with 

the question of robustness of the test to the unspecified 

parameters, namely, p and g , g . In order to determine 

X y 

the sensitivity of the test to the misspecified correla- 
tion coefficient p, the following distance was evaluated 
for different values of 

X y 

D^(p) = n l^^iGCe) - B(e)}^ dG( 0 ) (4.4) 

where G(e) is given in equation (3.2) when p = 0, and 
B(e) is given by equation (3.2) whenever p > 0, for 
p = .25 (.25). 75. Table 4 contains the type 1 errors 
given by 

Pr[W^ > D^(p)] = a^(p) (4.5) 

for different values of n and ~ ' “n^P^ 
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4 


Table 1 . Tail Probabilities for a ^( 6 )* 


r , 

n 

fi =.25 

.50 

.75 

1.00 

1.25 

1.50 

1.75 

2.00 

1 0 

20 

1.00 

.45 

.22 

. 12 

.06 

.04 

.02 

.02 


50 

.53 

.12 

.02 

.01 






100 

.25 

.02 

.01 






.25 

20 

1.00 

.41 

.18 

.09 

.05 

.03 

.02 

.01 


50 

.48 

.09 

.02 

.01 






100 

.21 

.01 







.50 

20 

.87 

.33 

.13 

. 96 

. 03 

.02 

.01 

.01 


50 

.41 

.06 

.01 







100 

.15 

.01 







.75 

20 

.69 

.21 

.07 

.02 

.01 

.01 




50 

.27 

.02 








100 

.06 








2 0 

20 

1.00 

.82 

.48 

.27 

.17 

.13 

.07 

.04 


50 

1.00 

.36 

.13 

.04 

.02 





100 

.60 

.13 

.02 

.01 

.01 




.25 

20 

1.00 

.72 

.40 

.22 

. 13 

.07 

.04 

.03 


50 

.91 

.29 

.09 

.01 






100 

.52 

.09 

.01 






.50 

20 

1.00 

.59 

.30 

.15 

.08 

.04 

.03 

.01 


50 

.73 

.21 

.05 

.02 






100 

.41 

.05 

.01 






.75 

20 

1.00 

.40 

.16 

.06 

-03 

.01 




50 

.52 

.09 

.01 







100 

.24 

.02 







3 0 

20 

1.00 1 

.00 

.66 

.43 

.28 

.18 

.12 

.08 


50 

1.00 

.55 

.25 

.11 

.04 

.02 

.01 



100 

.87 

.26 

.07 

.01 





.25 

20 

1.00 1 

.00 

.56 

.35 

.21 

.13 

.08 

.05 


50 

1.00 

. 46 

.18 

.07 

.02 

.01 




100 

.59 

.11 

.01 






.75 

20 

1.00 

.56 

.26 

.12 

.05 

.02 

.01 



50 

.72 

.18 

.04 

.01 






100 

.39 

.04 

.01 







*if < . 01 , then the entry is left blank . 
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(4.5) is the expected type I error as a function of the 
nuisance parameter p. It should be mentioned that the 


distribution for is not the same as that given by 

Anderson and Darling asymptotic approximation since the 

nuisance parameter p is unspecified [cp., Stephen (1976)], 

however, it does not appear feasible to follow Darling's 

procedure for computing the exact distribution whenever 

p and 6 , 3 are replaced by their consistent estimators, 

X y 

In spite of this shortcoming, equation (4.5) is used. 
However, Stephens (1976) showed that the asymptotic 
approximation given by Anderson and Darling is conserva- 
tive as compared to his fitted distribution in the family 
of normal distributions [Stephens (1976) Table 4, p..367] 
and the extreme value distributions [Stephens (1977) 


Table 1, p. 687]. Thus, it seems reasonable that 
equation (4.5) is also conservative, that is, if oIj^(p) 
is the true value for the l.h.s. of equation (4.5), then 
a(p) < 


Table 2. Type 1 Errors for Unspecified 


Y 

n 

p = .25 

.50 

.75 

1 

20 

1.00 

1.00 

.65 


50 

1.00 

.87 

.25 


100 

1.00 

.50 

.06 

2 

20 

1.00 

1.00 

.55 


50 

1.00 

.74 

. 18 


100 

1.00 

.41 

.04 

3 

20 

1.00 

1.00 

.53 


50 

1.00 

.70 

. 16 


100 

1.00 

.37 

.03 
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From Table 2, it follows that the procedure is only 
sensitive to p whenever p = . 75 and n > 50. This obsei^/a- 
tion was also supported in the simulation study. 

In order to determine the sensitivity of the test to 
the scale parameters, the distance given by 

Dn(s) = n/^''^ {G( 0 ) - C(e)}^ dG(e) (4.6) 

where G(e) is given in (3.2) and C(e) is given by (3.2) 

whenever tane = sr, s = = . 90 (. 02) 1 . 10 . Errors in 

X y 

either of the scale parameters can be considered by varying 
s in (3.2). Table 3 contains the expected type I errors 
given by 

Pr[W^ > D^(s)] = a^(s) (4.7) 

for different values of n, y , and pi 

From Table 3 one observes that the procedure appears 
to be resilient to errors in the scale parameter and that 
one might have a type I error when y = 3, p = .75, and n = 
100 at the 95% significance level. In addition it also 
appears that the results are symmetric about s = 1. 

5. CONCLUSIONS AND SUMMARY 

A procedure is outlined for determining whether 
a four or five parameter bivariate gamma model is appro- 
priate. The procedure was evaluated and three 
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Table 3 . 


Type 1 Errors for 
Misspecified Scales 



Y 


1 


2 


3 


p 

n 

.90 

.92 

. .94 

1.06 

1.08 

1.10 

0 

50 

1.0 

1.0 

1.0 

1.0 

1.0 

1,0 


100 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

.25 

50 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 


100 

.89 

1.0 

1.0 

1.0 

1.0 

.87 

.50 

50 

1.0 

1,0 

1.0 

1.0 

1.0 

1.0 


100 

.67 

.96 

1.0 

1.0 

1.0 

. 68 

.75 

50 

.69 

1.0 

1.0 

1.0 

1.0 

.73 


100 

.37 

.57 

.89 

.83 

.60 

.40 

0 

50 

.98 

1.0 

1.0 

1.0 

1.0 

1.0 


100 

.56 

.79 

1.0 

1.0 

.80 

.57 

.25 

50 

.78 

1.0 

1.0 

1.0 

1.0 

,80 


100 

.44 

.64 

1.0 

1.0 

.66 

.46 

.50 

50 

.58 

.81 

1.0 

1.0 

.84 

.61 


100 

.38 

.46 

.74 

. 76 

.49 

.31 

.75 

50 

.31 

.49 

.78 

.80 

.52 

.34 


100 

.09 

.22 

.44 

.46 

.24 

.11 

0 

50 

. 68 

.98 

1.0 

1.0 

1.0 

1.0 


100 

.36 

.55 

.86 

.88 

.57 

.38 

.25 

50 

.55 

.77 

1.0 

1.0 

.80 

.57 


100 

.26 

.43 

.70 

.72 

.45 

.28 

.50 

50 

.37 

.57 

.89 

.92 

.60 

.41 


100 

. 14 

.28 

.52 

.54 

.30 

.16 

.75 

50 

. 17 

.32 

.57 

.59 

.34 

.19 


100 

.03 

. 10 

.27 

.28 

.12 

.04 


The values of s = (. 96 , 1 . 04 ) V7ere omitted since the Type 1 
error was 1.0 for all parameters. Likewise^ whenever n = 20 . 
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different functions were evaluated in order to determine 
the procedure's feasibility and sensibility to the nuisance 
parameters. Admittedly, the evaluation is very limited 
and there are several limitations which would prohibit 
this type of procedure as an omnibus test of fit. How- 
ever, the results appear to be promising to the experi- 
mentalist interested in obtaining insight into the stated 
problem. 

There are several nonparametric procedures for test- 
ing (’1.2) versus (1.3) and perhaps these are not as sensi- 
tive to the nuisance parameters. However, the proposed 
procedure is based upon "measuring" significant departures 
of the parametric distributions function which are vital 
to the modelers' primary objective. 
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CHAPTER III 


A DIFFERENTIAL EQUATIONS APPROACH TO THE 
MODAL LOCATION FOR A FAMILY OF 
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ABSTRACT 

Analytical and niomerical computational methods are 
given for determining the location of the mode as a func- 
tion of the parameters of a class of the bivariate gamma 
distribution. 


I . INTRODUCTION 

Smith, Adelfang, and Tubbs (1983) derived some compu- 
tational results for a family of bivariate distributions. 
In their paper they consider the location of the mode as a 
function of the shape parameters, and y 2» ^nd the 
dependence coefficient n . The purpose of this paper is 
to consider this problem in greater detail. Tha* is, the 
paper will consider analytical and numerical computational 



methodt. for locating the modal values for the class of density 
functions given in Smith, Adelfang, and Tubbs (1983). The 
general density function is given by: 


f ( i » ^2 * 


where a 


jk 


-(Cj^+t2> 

ti t£ e 

1.Y2.H) = 

(1-n) r (Y]^)r (y 2 -Y]^) 

(Y2"Yi+k) 

(l-n)^^‘*‘^r(Y2+j+k)j:k; ’ 


z z 

j=»0 k-0 



( 1 . 1 ) 


tf = Bj^x, t 2 = S 2 y. &i,&2 scale parameters, Y 2 > Ti > 1 

are shape parameters, and 0 < n < 1 is associated with the 
correlation coefficient p by the equation n = p^y'^TT^. We 
will assume without loss of generality that ~ ^2 ~ 

We will concentrate on the special case = Y 2 = Y of 
(1.1) for which the distribution function reduces to 


f(t^,t2iY f^) 




_l - (t^+t 2 )/ (1-n) 


(tj^t2)^ 


(l-n)T(Y) 


j=0 (l-n)^^r(Y+j)j: 

( 1 . 2 ) 


This is the form given by Kibble (1941) . 


Smith and Adelfang (1981) used the above class of density 
functions in modeling wind gust data for the ascent flight of 
the Space Shuttle. A parametric model was selected in that the 
parameters are used to establish engineering constraints for 
the shuttle payload system. Thus, the modal location and value 
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were of interest to this particular application. The authors 
are not aware of any other results, either analytical or numer- 
ical, for the modal location for non-Gaussion multivariate dis- 
tributions. The closest related work is in the area of density 
and mode estimation [3.g. Sager (1978, 1979), de Beauville 
(1978), and Eddy (1980)], 

In Section 2 we will derive some qualitative results con- 
cerning the behavior of the modal location of (1.2) as a func- 
tion of (y.h). Section 3 presents analogous results for 
another borderline case = 1» Y 2 = ^ of (1.1). In Section 
4 we present some numerical procedures based on the theoretical 
investigations of the previous sections. The general case 
Y 2 > Yi > 1 is considered in Section 5. We present some numer- 
ical tabulations for the modal location of (1.1) as a function 
of (y]^, Y 2 » consider some numerical interpolations from 

■the borderline cases considered in Sections 2 and 3. 


2. EQUAL SHAPE PARAMETERS - ANALYTICAL METHODS 

Lemma 1 . The function f (t^^ , t 2 ; y . n) defined by (1.2) attains 
its maximum in the region R;^ = { (t , t 2 ) : t^^O , t 2 ^ 0 } on the 
line t^ = t 2 * 

2 

Proof; Since f is integrable and continuous over R_^, it is 

2 

clear that f attains its maximum on R_j_. Choose any constant 
c > 0. Let h(t) = f (t,c-t;Y ,n) , 0 < t < c. Then from (1.2) 
we have 
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h(t) - I 

j=0 J 

where Kj(Y,n) > 0 is independent of t. Therefore, 

h'(t) = r K. (Y,n)(Y+j-l)t^'^-"^(c-t)^'^j‘2(c-2t). 
j=0 

Since h'(t) > 0 for 0 < t < c/2 and h'(t) < 0 for c/2 < t < c, 
h attains its maximum at t = c/2.. Therefore f(t^,t 2 ) attains 
its maximum along any line tj^+t 2 =c at the point (c/2, c/2). 

This completes the proof. 

Define g(t;Y,n) » f(t,t;Y,n). Then by Lemma 1 it is suf- 
ficient to find the point on t ^ 0 at which g attains its maxi- 
mum value. Using (1.2) one can show chat 

g(t;i-,n) - (2.1) 

where c(Y,n) = [ (1-n) (*^)'^"^r(Y)J and h(t) = ^ .(p(ri)c). 

Y 

where I (z) denotes the modified Bessel function with index u , 

V 

and p(n) = 2/n”/(l-n). 

Using [Abramowitz and Stegun (1964), Eqn. 9-6-28] it is not 
difficult to show that h' (t) = p(n)t^ i,y_2 ' therefore 

f (.j ^ j,) is the mode at the bivariate gamma distribution given 

by (1.2) if and only if g'(t) = 2g(x)/ (1-n) or 

*^^Y-2 = ly-i^P ('i)t) . (2.2) 

where p(n) = 2/r7/(l-n)- 

With the aid of (2,1), we may prove the following theorem. 

Theorem 1 . For fixed y > 1, let t(ti) denote the value at which 
f(T(n), T(n);Y.n) is a maximum. Then t is continuously differ- 
entiable for 0 ^ Tj i ^’od satisfies the initial value problem 
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x'(n) - (T/2n)((2T- 2^+3)"^ - •+n)(l-n)’b 

(2.3) 

t(0) = Y-1. 

Proof; It is easy to show directly from (2.1) that g attains 

its maximum at t = Y - 1 when n. = 0# so that t (0) = y ~ 1* 

Furthermore is continuously differentiable for 0 < n < 1 and 

computation shows that ~ — 2. ^ 0 at t = yi and n=0. Therefore, 

3t2 

T(n) is continuously differentiable in a neighborhood of n = 0 
by the implicit function theorem. The proof will be completed 
by differentiating both sides of (2.1) with respect to n • 

After some simplification and solving for T'(n) this yields 

= (l-n)l 5 (p (n) T )/ (4nq (n ) ) - (1+n ) x/2n (1-n ) (2.4) 

where q(n) = I^-l ^ ^ 1^-2 (p ^ ^ • 

By [Abramowitz and Stegun (ly64) , Eqn. 9-2-26], 

^Y-i^P “ (Y“i)I^_3^(p(n)x)/p(n)T 

= ^Y_i(P(’i)t) + (y-2) I ^_2 (p(n) x)/p(n)x ; 

substituting these expressions into q(n) and using (2.2) .yields 
after some simplification 

q(n) = (1-n) (2t-2y + 3)I^__2 (P(’^Jt)/2t. 

Substituting this expression into (2.4) completes the proof 
of Theorem 1. 

The nonlinear differential equation (2.3) cannot be solved 
in general in closed form. Some numerical solutions arc given 
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in Section 4. However (2.3) does give information regarding 
the qualitative and limiting behavior of T(n) for y > 1, In 
the special case y = 3/2, (2.3) reduces to a linear differen- 
tial equation which can be solved directly by standard methods. 

Corollary 1 If y = 3/2, then 



This result can also be obtained directly from (2.2) using the 
fact that (z) can be expressed in terms of hyperbolic func- 
tions when u = ±1/2. 

Since the differential equation (2.3) is singular at 
n = 0, its numerical solution requires some additional knowl- 
edge of the behavior of x(n) near n = 0. This is provided 
by the following corollary. 

Corollary 2 . The function t(ii) is continuously differentiable 
at n = 0 and satisfies 

t’(0) = -(y-l)/Y . Y > 1. (2.5) 

Proo f. The continuous differentiability of t at n = C was 
considered in the proof of Theorem 1. Choose n > 0 , then by 
the mean value theorem there is a number 5 e (0,n) such that 
±(n) = t(0) + ht’( 5) = y-1 + nx'(5). 

Substituting this expression into (2,3) and simplifying yields 

T ' (n) = (Y-l-fm ' (K)) (l+(l+n)T * (C)) 

(l-n)(l+2riT '(O) 
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Letting n 0 and using the continuity of r'Cn) we have 








T'(0) = -(y-1)(1+t’(0)). 

Solving this equation for t ' ( o) yields (2,5) 

We will write when we wish to emphasize the aepend- 

ence of the modal location on y* Theorem 1 and Corollary 2 
may be used to obtain several of the qualitative and asymptotic 
properties of the function T(n»Y) the region 0 < n < 1/ 

Y > 1. These are summarized in the following theorem. 


Theorem 2. The modal location function has the following 

properties : 

(i) x(n,Y) is a decreasing function of n for fixed y > !»* 

3 

(ii) lim T(n»Y) = max (y-^.O- fot y > 1; 

n^l ^ 

3 

(iii) T(nfY) ~ (y~j) is a decreasing function of y 
fixed n e (0,1) and y > 1; 

(iv) lim T(n,Y) -(y“^)= 2Ti+^* 0 4 n 4 1. 


Proof : We will show that x' (n) < 0 for 0 4 n < 1- Suppose 

not, then since x * (0) < 0 by Corollary 2, there is a point 
5 > 0 such that t ' (5 ) = 0 and x ' (n) 0 for -0 < n < C • Let 

w(n) = x(n) - (y-|-) and z(n) = 2lTTnT' from (2.3) it is 

easy to see that 


X ' (n) 


X (n) r 1 _ 

4n ^w (n ) z (n ) “^ ' 


so X ' (5 ) = 0 if and only if w (£ ) = z { 5 ) . 
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Let h = w - z. Note that z'(n) = -(1+n)’^ so that h'(0) = 

w'(0) - z' (0) = t' (0) + 1 > 0 and h(0) = w(0) - z(0) =» 0. 

Therefore since h(C) = 0 and h(n) > 0 for 0 < n < £ , we must 

have h'(5) 4 0. However, h'(5) * w'(t) - z' CO = t'(C) - 

z'(C) = -z'(C) = (1+C)”^ > 0. This contradiction proves (i) . 

Furthermore, we have that w(n) > z(n) for 0 < n < 1. 

We will now consider the proof of (iii) . Fix Yi Y 2 ^ ^ 

and let f(n) = w(n,Y]^) - w(n,Y2) before w(n,Y) = 

■c(n,Y) - (y“^)- We wish to show that f(n) < 0 for 0 < n < 1. 

Clearly f(0) = 0 and by (2.5) f(0) = i i <0. Assume 

^1 ^2 

to obtain a contradiction that there is a point Ce( 0,1) such 
that f(5) = 0. If, in addition, we assume C is the first 
such point, then f(n) < 0 for 0 < n < 5 so f '(c) i 0. However, 

using (2.6) at both and Y 2 the fact that w(C,Y]^) = 

w(c,y 2 ) ft is not difficult to show that 

(yi-Y2> 


f'(0 = 


'+T1 


1 1 
w(c) ' z(c) 


Since y^ > Y 2 ^ z(5), it follows that f'(c) < 0. 

This contradiction completes the proof of (iii) . 

Now we turn to the proof of (ii). First consider the 
case 1 < Y 4 3/2. Since t is decreasing in n and positive 
for 0 < n < 1 we know that t* = lim t(ti) exists, where the 

“ T-*-i 

limits at 1 are always from the left. Assume to obtain a 
contradiction that t* > 0. Then it is not difficult to show 
using (2.3) thit 
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w ^ 1 T*(l+n) 

" ^ ■ 2rr(I-TT 

Therefore, for j ^ n < 1 we have 
= i ■ 2 (l-n) • 

Integrating both sides of this inequality from ^ to n yields 
t(ti) 4 T (^) + j ^ In(l-n) 

for ^ 4 n < 1. However, this implies that x(n) -► - » as 
n ■+■ 1 , a contradiction. 

The case y > j follows easily from (iii) and the proof of 

3 

(i) because for y ^ 2 " 

z(n) <, x(n) - (y-^) 4 

O 

and both z(n) and xCn.j-) approach zero as n -*• 1 . 

Finally, we consider the proof of (iv) . Let u(ti,y) “ 

3 

w(n,Y> - z(n) for 0 < n < 1 and y > Fi’oid the proof of (i) 

we know that u(ti,y) ^ 0- From (2.o) we obtain 

TT» ^ ^ n - ^(n .y) 1 . r 2 (t|) - w(n ,Y) 

w (h.y) 5tT [_ ~ z (n) J *"5n L w(n,Y)z(n) _ 

so that 

w'(ti,y) < - (Y-f)vi(n ,y) • 

Therefore, 

u'(n,Y) “ w'(n,Y) - z'(n) 4 -(y-|)u(ti,y) + 1. 
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0 



From this inequality we obtain 




u(n - u( 0 ,y) 4 — -1 . 

(y - j ) ^ -* 

Therefore 

0 < u(ti,y) 4 — 

This implies that u(n,Y) 0 as y “ and compleces the 
proof of Theorem 2. 


3. UNEQUAL SHAPE PARAMETERS- -THE CASE Y^ = 1 


In this section we consider another "borderline” case of 
the general bivariate gamma distribution, the case Yj^ * 

For technical reasons we will limit our discussion to the 
range Y 2 i 2 and for brevity let Y 2 “ y- Then the function 
given by (1.1) reduces to 


f(tj^,t2;l,Y,n) 


si 

(l-n)f(Y-lT jfg® ^ F 


(3.1) 


where 


c. 

J 


£ 

k=0 



r(Y-Hc-l) . 

kir (Y+j+k) ’ ^ 


0,1,2... 


ti t2 

and where si = t — , s o “ t — . and s, 
i 1-n 2 i-n 3 


ns2. 


(3.2) 


The following lemma allows us to restrict our attention 
to the line tj^ = 0. 
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0 



Lemma 2 . The function f (tj^, t 2 ;l,Y. n) given by (3.1) for 
Y i 2, takes on its maximtmi value in the region ^ 0, 
t 2 i 0 on the line t^ = 0. 

Proof: Since f is continuous and integrable in the first 

quadrant, we know it takes on its maximxim value at some point 

“(fc 

(t 2 ^,t 2 ). We will prove that tj^ = 0 by showing that for any 
fixed t 2 > 0 , f(tj^,t 2 ) is a decreasing function of This 

is equivalent to showing that the function 

” j 

g(s) = z e“® c. Iy 
j=0 J 


is a decreasing function on a ^ 0 where c. is given by (3.2). 
Note that 


g* (s) = -e"^ Z c. + e"® Z c, 


L c. T e i c. 
j— Q J J j=0 ^ 


. -e-5 i c. 4 + e'® £ =•+1 4 

j.O J F J.0 Jl 


. .e-s ~z 4 


4oF 


Therefore g'(s) < 0 for s ^ 0 if < c^ for j = 0,1,2 

To this end note that 


c 


j+1 


z 


k=0 


j+l+k r(y-Kk-l) 

®3 Rlr (Y+j+l+k) 


00 

z 


k=l 



r(y+k-2) 
(k-1) !r(Y+j+k) 


= E r(Y+k-l) . . k < 

®3 klr (Y+j+k) PTn 


c. 

J 
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since y i 2 implies chat <, 1 for k » 1,2,3,.... This 

completes the proof of Lemma 2. 

According to the preceding lemma, the mode of the bivariate 
gamma distribution in this case is the point (0,jj) where m is 
the point on t ^ 0 where the following function is a maximum: 


g(t) » t^"^ h(t) 


where 


h(t) 


/n t s r (v+k-l) 


kio kl(-r+tc-l) 


Note that 


h(t)) 


“ ^ Y"^h-1 


--n 


k»0 


-n^ kl (y+k-l) 


) 


_ n • ,nt ^ 

- tttt <Trrr> kr 


= ^ gnt/(l-n) 


Therefore, 
h(t) = 




■n o 'i-ir 




SO that the function we wish to maximize is 


g(t) = s^*'^ gHs/Cl-n) ds, 


(3.3) 
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Lennna 4 . Let y(n), or when necessary y(n,Y). denote the value 
for which f(0,y(n,Y);l,Y.n) is a maximum where f is defined 
by (3.1) and (3.2). Then y (0) = Y“1. v (n) i Y"2 for 0 4 n < 1» 
and satisfies the eqiiation 

g(y) = (l-n)y^ ^ , 0 ^ n < 1, (3.4) 


where g is defined by (3.3). 

Proof ; It is easy to see from (3.3) that g attains its maximum 
on 10,-) at a point t* > 0 for which g'(t*) = 0 and g"(t*) ^ 0. 
Differentiating (3.3) we obtain 


g'(t) = - g(t) + t^'^e"'^ 

and 

g"(t) - - g'(t) + 


Therefore g' (y) = 0 implies (3.4) and g''(y) ^ 0 implies that 

•t t^ 


y Y " 2. Since when 


n 


0, g(t) = e 


•y-T 


It is easy to 


see that y (0) = y~1. This completes the proof of the lemma. 

With the aid of these preliminaries we may prove the 
following theorem in the spirit of Theorem 1. 


Theorem 3 . For fixed y i 2, let y(n) denote the value of 
which f(0,y(n); l.Yin) is a maximum. Then y is continuously 
differentiable on 0 ^ n < 1, y ' (0) = -1 + ^, and on 0 < n < 1 
y satisfies the initial value problem 
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(3.5) 


y(0) - Y-1. 

Proof ; As in the proof of Theorem 1 the continuous differen- 
tiability of jj in a neighboi'hood of n=0 may be proved by applying tht 
implicit function theorem to (3.4). This differentiability will 
be extended to all of [0,1) by proving that (3.5) holds. Let 

g(t,n) denote the function defined by (3.3) and let g and g 

t n 

denote its partial derivatives with respect to t and n , 
respectively. Then differentiating both sides of (3.4) with 
respect to n we obtain 

g|-(iJ,n)y' + g^(y.n) = (l-n)e ^ (y-2-jj )ji ' -e (3.6) 

By definition g^(y,n) * 0 and direct differentiation of (3.4) 
and integration by parts yields for 0 < n < 1 that 

+ e-c/a-n) jY-l eis/(l-n)j3 

/i _\2 O 


j g(t,n) 

(1-n)^ 


+ 


1 

7 

(1-n)^ 


-t/(l-ri;' r 

® ^ n 


^Y-lgnt/ (1-n) 


- (y-1) /q s^‘V®/^^"'’^ds] 
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0 



(l-n 








Therefore, using (3.4) we obtain 


g^(w.n) 


U Y“2 -U , 1 Y-1 “li 

c— y ' ‘e ^ + — >4 V V ^ 

-n n(l-n) 


•1 y~ 2 -y 

— y ' e 


. + >> . izi) 

l-n n(l-n) n ' 


we ^ . 

Substituting this expression into (3.6) and simplifying yields 
(3.5). For n =* 0, an easy calculation shows that 

from which substitution into (3.6) with n “ 0 and y * y “ 1 
shows that y ' (0) = -1 + ^. This completes the proof of 
Theorem 3. 

The following corollary exploits the fact that (3.5) 
reduces to a linear differential equation when y “ 2. 

Corollary 3 . If y * 2, then 

y(n) * ^ 

Proof . For y • 2 equation (3.5) reduces to 

. . y _ 1 

^ rTn--n) n 
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I* 

f' 



which is easily solved in closed form by standard methods to 
show the desired result. This result is also easily derived 
directly from (3.4). 

It is interesting to note that the translated modal 
location function v(n) ■ y(n) - (y"2) satisfies the differ- 
ential eqioation 


v’(n) 


v(n) + y-2 
n 


( 


1 1 V 


v(0) 


1, v’(0) 



whereas the translated modal location function w(o) =* 

3 

T(n) - of Section 2 satisfies the analogous differential 

equation 


w' (n) 


w(n) + Y 


3 

7 




( 1 2(l+n)x 


w(0) - ^ , w' (0) - -1 + 

For this reason y behaves in a manner similar to t. Its 
properties are stated in the following theorem. Since the 
proof of this theorem is entirely analogous to tne proof Yf 
Theorem 2, it is omitted. 


Theorem 4 . The modal location function y(n,Y) the follow- 
ing properties: 

(i) y(n,Y) is a decreasing function of n for fixed 
Y i 2; 

(ii) limy(n,Y) " y " 2 for y i 2; 
n-*-l 
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0 



(iii) y(n,f) - (y~2) is a decreasing function of y for 
Y ^ 2 and fixed n e (0,1); 

(iv) lim(y(n,y) - (y" 2)) = 1-n for 0 4 n 4 1. 

Y ->-00 

4. NUMERICAL RESULTS 

In this section we present some quantitative results based 

on the results of the previous sections. Table 1 shows the 

value of the modal location function for equal shape parameters 

T(ntY) various values of n and y. Table 2 shows values of 

3 

the translated modal location function w(ti,y) = tCh.y) ~ (y“^)* 
This table illustrates the qualitative behavior of this fianc- 
tion derived in Theorem 2. The limiting values of n = 1 and 
Y = “ are taken from Theorem 2 . 

The values in Tables 1 and 2 were computed using Theorem 
1, Specifically, a fourth-order Runge-Kutta algorithm was used 
to compute an approximate solution of the differential equation 
(2.3) on the interval 0 4 n < 1 for each specified value of y. 
Since equation (2.3) is singular at n = 0, Corollary 2 was used 
to replace the initial condition t ( 0 ) = Y~lhy the approximate 
initial condition 

T(h) » T-1 - h 

where h is the step size of the numerical method. Figure 1 
shows the data of Table 2 in graphical form and illustrates the 
behavior of the function w(n,y) derived in Theorem 2. 
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Tables 3 and 4 show the corresponding results for the 
modal location function y(n,Y) the case 

and its associated translate v(n,Y) = u(n.Y) " (y"2). These 
tables were computed by the same methods as Tables 1 and 2 
except using the results of Section 3. Figure 2 illustrates 
the qualitative behavior of the function v(n,Y) as indicated 
by Theorem. 4. 

Note that the differential equations (2.3) and (3.5) 
allow the modal location to be computed recursively in n for 
a fixed value of y as a dynamic process in a time scale meas- 
ured by the modified correlation coefficient n . Error in the 
computation is introduced through the discretization of this 
continuous evolutionary process. A more conventional compu- 
tation of the modal location would require an independent 
calculation for each value of n with error introduced through 
the truncation of the series representation (1.2) of the dis- 
tribution function. This error becomes particularly trouble- 
some as n 1. 
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Tab? j 1. Selected values of the modal location function T(n.v) 
for equal shape parameters. 


Y\n 

0 

. 1 

.3 

.5 

.7 

.9 

1 

1. 1 

.1000 

.0908 

.0720 

.0525 

.0322 

.0110 

.0000 

1.3 

.3000 

.2765 

.2268 

.1724 

.1117 

.0412 

.0000 

1.5 

.5000 

.4660 

.3931 

.3116 

.2169 

.0958 

.0000 

2.0 

1.0000 

.9491 

.8412 

.7277 

.6127 

.5279 

.5000 

3.0 

2.0000 

1.9334 

1.8034 

1.6862 

1.5939 

1.5268 

1.5000 

10,0 

9.0000 

8.9152 

8.7754 

8.6697 

8.5891 

8.5624 

8.5000 


Table 2. 

Selected values of the 
w(n,Y) for equal shape 

translated modal 
parameters . 

location 

function 

y\n 

0 

.1 

. 3 

. 5 

.7 

.9 

1 

1. 1 

.5000 

.4908 

.4720 

.4525 

.4322 

.4110 

.4000 

1.3 

.5000 

.4765 

.4268 

.3724 

.3117 

.2412 

.2000 

1.5 

.5000 

.4660 

.3931 

.3116 

.2169 

.0958 

.0000 

2.0 

.5000 

.4491 

.3412 

.2277 

.1127 

.0279 

.0000 

3.0 

.5000 

.4334 

.3034 

.1862 

.0939 

.0268 

.0000 

10.0 

.5000 

.4152 

.2754 

.1697 

.0891 

. 0264 

.0000 

00 

,5000 

.4091 

.2692 

. 1667 

.0882 

.0263 

.0000 
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5. UNEQUAL SHAPE PARAMETERS 


In this section we briefly consider the mode of the 
general bivariate gamma distribution given by (1.1). By 
setting the partial derivatives of f (t-j^, t2;Yi»Y2’'^^ with 
respect to tj^ and t£ equal to zero, one finds that f attains 
its maximum at the point (tj^,t 2 ) whose coordinates satisfy 


t. 



I 

j=0 




(5.1) 


and 


to 




Z Z a.. (Y 2 +j+k-l) 
c=0 j=0 


(5.2) 


00 

z a.. 


where S = Z ^ 
k=0 j=0 

and a^j^ given as in (1.1) depends on t^^ and t 2 * 

Table 5 shows selected values of the modal location for 
the case Y£ “ 3. They were computed by truncating each of 
the series in (5.1) and (5.2) to about fifty terms and simul- 
taneously iterating on these equations until an approximate 
solution is obtained. These computations become unreliable 
as n 1 and the truncation error becomes tonaccep table. 

Figure 3 gives a graphical representation of the change 
in modal location with n and for fixed y£ “ It is 

interesting to note for a fixed n the extent to which the 
modal location may be approximated by linear interpolation 
between the borderline cases discussed in Sections 2 and 3. 


III-23 


More specifically, vre have the empirical approximations 




^1 


Yi-1 


(5.3) 


and 


t 


2 


V (n , Y£) + 



(t (n ,Y2) 


y (n ,Y2)) 


(5.4) 


where t and y are as defined in Sections 2 and 3, respectively. 
This empirical relationship is a subject for further investiga- 
tion. 
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CHAPTER IV 


ANALYSIS OF WIND GUST DATA 
J. D. Tubbs 

Department of Mathematical Sciences 
University of Arkansas 
Fayetteville , Arkansas 

ABSTRACT 

This paper summarizes the analysis of wind gust data 
using statistical and mathematical procedures which were 
developed for the bivariate gamma distribution. 

1 . INTRODUCTION 

Adelfang and Smith (1981) discuss the use of the gamma 
distribution in modeling gust data at Cape Canaveral, 
Florida. Smith and Adelfang (1981) treated gust amplitude 
and length scale as the variables of the bivariate gamma 
distribution. Smith, Adelfang, and Tubbs (1983) presented 
some useful analytical and computational results for a class 
of the bivariate gamma and applied some of these results to 
the wind gust data. The purpose of this paper is to analyze 
the wind gust data using some additional analytical results 
obtained for the bivariate gamma distribution. 

2 . DATA 


The data used in this paper consists of absolute gxist 


magnitude and gust length for both the zonal and meridional 
components. Tha 130 wind profiles were filtered using the 
band pass filter for wavelengths within 420-2470 meter band. 
Data were available for the reference altitudes: 4Km, 6Km, 
8Km, lOKm, 12Km, and 14Km. The data was paired into bivar- 
iate components for both the zonal and meridional components, 
denoted by the pairs (Au,Lu) and Av,Lv), respectively. 

3. ANALYTICAL PROCEDURES 

The data were partitioned according to reference alti- 
tudes, then the 150 observations were analyzed using both 
univariate and multivariate techniques. Simple descriptive 
univariate techniques were generated using PROC UNIVARIATE 
in SAS. These procedures were used to help in the assess- 
ment of the marginal distribution. The multivariate 
descriptive procedures consisted of bivariate scatter plots 
and contour plots. 

Goodness of fit tests consisted of a univariate test 
for marginal nozmality generated by SAS, two tes^s for 
bivariate normality as discussed in Meredith and Tubbs 
(1981), and a bivariate test for the gamma distribution. 

The latter procedure is a bivariate Chi-square type test 
which uses the computational results for the distribution 
function as presented in Smith, Adelfang, and Tubbs (1983). 

Parameter estimates for the bivariate gamma distribu- 
tion were evaluated. These estimates were then used in 
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generating the three-dimensional bivariate gamma density 
function plots and the modal locations were estimated 
using the results given by Brewer, Tubbs, and Smith (1983). 

4. RESULTS 

The results are summarized in Tables 1-7. Additional 
results are given in Appendices A and B. 

Tables 1-6 summarize the results for both the test of 
fit and the parameter estimates for the bivariate gamma 
distribution. There are two main tests for bivariate nor- 
mality and both of these are discussed in Meredith and 
Tubbs (1981). The first is a procedure proposed by Rincon- 
Gallardo et al. (1979). Since this procedure transforms 
the data to a univariate test for uniformity three differ- 
ent tests for uniformity are used. The second test for 
bivariate normality is based upon a procedure proposed by 
Cox and Small (1978). 

The bivariate test for the gamma distribution is a Chi- 
sqviare type test of fit. Thus, this procedure has the 
usual difficulties of selecting the number of cells and 
cell location that are associated with this type of test. 

In the interest of time and space a fixed procedure was 
applied for all the data sets. Namely, the marginal distri- 
butions were partitioned according to the .05, .25, .50, .75, 
and .90 quantiles based upon the gamma parameter estimates. 
This particular choice affected the results for some of the 
data sets, however, it seemed a reasonable global choice. 
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The univariate gamma parameters were estimated using 
a maximvim likelihood procedure presented by Greenwood and 
Durand (1960) and discussed in Tubbs and Brewer (1981). 

Appendix A contains the results for the univariate 
descriptive statistics. Appendix B contains plots for each 
data set. The density functions were generated using the 
gamma parameter estimates. The contour plots are level 
slices of the density fxinction and are not equal proba- 
bility contours. The location of the mode is denoted by 
the symbol + and this value is computed using the analytical 
results given in Brewer et al. (1983). Table 7 summarizes 
the results for the v /dal location. 
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Table 1. Sunmary for Wind Gust Statistic 

Using Band Filter 420-2470 Altitude - 4 Km. 


Multivariate Test 



(Au , Lu) 

(Av.Lv) 

Cramer-Von 

Mises 

.2062 

.2144 

Normality Watson's 


.2023* 

.2058* 

K - 

S 


.0618 

.0551 

Cox 



11.45** 

26.07*** 

Gamma Chi- 

■Square 


33.3 

53.00*** 

Univariate Test 





Au 



.077* 


Normality 



.068 
. 090** 


Lv 



.077* 


Parameter Estimates 


A 

Y 

e 

p 

Au 


3.402 

2.430 

.2280 

Lu 


5.275 

.007 


Av 


2.808 

1.891 

.3415 

Lv 


4.429 

.006 



* denotes 

** denotes 

*** denotes 


that test is 
that test is 
that test is 


significant at 
significant at 
significant at 


.05 level. 
.01 level. 
.001 level. 
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Table 2. 


Surmnary for Wind Gust Statistic 

Using Band Filter 420-2470 Altitude « 6 Km. 


Multivariate Test 



(AuLu) 

(Av.Lv) 


Cramer-Von : 

Mises 

.3806 

.2942 

Normality 

Watson's U^ 


.2411* 

.2208* 

K - 

S 


.0897 

.0623 


Cox 



24.4*** 

31.8*** 

Gamma 

Chi- 

Square 


72 . 03*** 

57.14*** 

Univariate Test 

Au 



.047 


Normality 

Lu 

Av 



.081* 

.054 



Lv 



.072 


Parameter Estimates 


A 

Y 

A 

0 

A 

P 


Au 


2.577 

1.916 

.4217 


Lu 


4.102 

.005 



Av 


3.168 

2.030 

.2506 


Lv 


4.895 

.005 



* denotes 

** denotes 

*** denotes 


that test is 
that test is 
that test is 


significant at 
significant at 
significant at 


.05 level. 
.01 level. 
.001 level. 
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Table 3. 


Sxunmary for Wind Gust Statistic 
Using Band Filter 420-2470 Altitude 


8 Km. 


Multivariate Test 



(Au , Lu) 

(Av.Lv) 


Cramer-Von 

Mises 

1.090*** 

.490* 

Normality 

Watson's U2 


. 721*** 

^ 409*** 

K - 

S 


.102* 

.083 


Cox 



10.67** 

8.67* 

Gamma 

uni- 

•Square 


34.74 

52.57** 

Univariate Test 

Au 



.114** 


Normality 

Lu 

Av 



. 108** 
. 107** 



Lv 





Parameter Estimate j 


Y 

6 

p 


Au 


3.023 

2.113 

.3396 


Lu 


4.149 

.005 



Av 


2.922 

1.811 

.4253 


Lv 


4.61«h 

.005 



* denotes that test is significant at .05 level. 

** denotes that test is significant at .01 level. 

*** denotes that test is significant at .001 level. 




Table 4. Summary for Wind Gust Statistic 

Using Band Filter 420-2470 Altitude = 10 Km. 


Multivariate Test (Au.Lu) (Av,Lv) 


Cramer- Von Mis es 

.129 

. 753** 

Normality Watson's 

.126 

.469** 

K - S 

.052 

.104* 

Cox 

8.23* 

8.58* 

Gamma Chi-Square 

45 . 04** 

36.29 

Univariate Test 



Au 

.073* 


Normality 

.063 

.073* 


Lv 

. 109** 


Parameter Estimates y 

B 

p 

Au 3.04]' 

1.909 

.4345 

Lu 5.203 

.006 


Av 2.522 

1.300 

.3890 

Lv 4.543 

.005 


* denotes that test is significant 

at .05 level. 


** denotes that test is significant 

at .01 level. 


*** denotes that test is significant 

at .001 level. 
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Table 5. Summary for Wind Gust Statistic 

Using Band Filter 420-2470 Altitude = 12 Km. 


Multivariate Test (Au.Lu) (Av.Lv) 



Cramer-Von Mises 

.465* 

^ nr 

.398 

Normality 

Watson's 

.391** 

. 329** 


K - S 

.077 

.076 


Cox 

25.04*** 

23.70*** 

Gamma 

Chi-Square 

52 . 54** 

41.24* 


Univariate Test 




Au 

. 116** 

Normality 

Lu 

Av 

. 112** 
.042** 


Lv 

. 085** 


Parameter Estimates 

Y 

A 

6 

P 

Au 

2.155 

.969 

.2983 

Lu 

4.113 

.005 


Av 

2.612' 

1.051 

.3066 

Lv 

4.023 

.005 



* denotes 
** denotes 
*** denotes 


that test is 
that test is 
that test is 


significant at 
significant at 
significant at 


. 05 level . 
.01 level. 
.001 level. 
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Table 6. Summary for Wind Gust Statistic 

Using Band Filter 420-2470 Altitude * 14 Km. 


Multivariate Test 


Normality 

Cramer-Von Mises 
Watson's u2 
K - S 

1.359*** 
. . 647*** 
. 134** 

.616* 

.357** 

.084 


Cox 

14.6*** 

17 _ 9 *** 

Gamma 

Chi-Square 48 

48.29** 

42.97* 


Univariate Test 




Au 

. 090** 

Normality 

Lu 

Av 

.096** 
. 085** 


Lv 

.061 


Parameter Estimates 

A 

Y 

e 

A 

P 

Au 

3.803 

1.549 

.2420 

Lu 

4.725 

.005 


Av 

3.324 

1.161 

.3171 

Lv 

5.057 

.006 



* denotes that test is significant at .05 level. 

** denotes that test is significant at .01 level. 

*** denotes that test is significant at .001 level. 
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Table 7. Modal Location 


Altitude 

Variables 


Method’^*’ 





I 


II 

4000 

(Au,Lu) 

.939 

581.9 

.939 

581.7 


(Av,Lv) 

.874 

516.7 

.873 

516.7 

6000 

(Au,Lu) 

.747 

539.6 

.745 

539.4 


(Av,Lv) 

1.008 

729.4 

1.007 

729.4 

8000 

(Au , Lu) 

.875 

571.0 

.874 

570.8 


(Av,Lv) 

.965 

639.8 

.960 

639.8 

10000 

(Au,Lu) 

.984 

624.8 

.977 

625.2 


(Av,Lv) 

1.064 

624.8 

1.058 

624.8 

12000 

(Au,Lu) 

1.083 

556.8 

1.080 

556.8 


(Av,Lv) 

1.402 

546.2 

1.400 

546.2 

14000 

(Au,Lu) 

1.713 

705.4 

1.713 

705.4 


(Av,Lv) 

1.868 

626.0 

1.865 

626.0 

* Method I 

Truncation of a 

double i 

series. 



Method II 

Interpolation. 







5. SUMMARY 


The data sets are discussed according to reference 

altitude. 

4 Km. The normal or the gamma are not rejected for the 
zonal (u) components. As discussed in Meredith and 
Tubbs (1981) the Cox and Small procedure is sensitive 
to symmetry and is not recommended for this data. The 
gamma distribution was rejected for the v-component and 
normality was accepted. However, marginal normality 
was rejected at the .01 level for the absolute gust 
magnitude (Av) . 

6 Km. The bivariate gamma was rejected in both wind com- 
ponents and normality was not rejected. 

8 Km. Normality was rejected for both wind components. The 
bivariate gamma was accepted in the u-component but not 
for the v-component. 

10 Km. The u-component appears to be normal whereas the 
gamma is accepted in the v-component. 

12 Km. Both distributions appear to be suspect for the 

u-component and the gamma is accepted for v. Normality 
is also rejected for v by considering the marginal dis- 
tributions. 

14 Km. Neither distribution is acceptable for u and the 
gamma is perhaps better for v. 
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6. CONCLUSIONS 


The wind gust data was analyzed using some new proce- 
dures for the bivariate gamma distribution. The analysis 
was meant to be informative, in that it represents examples 
for some of the analytical procedure. The analysis is not 
meant to be completel * *-horough. Hence, there are still 
some unresolved questions concerning the applicability of 
the bivariate gamma for modeling wind gust data. One sus- 
pects that neither the normal nor the gamma are completely 
appropriate, however, perhaps both could provide acceptable 
results for defining engineering constraints. 

As mentioned in the paper the test for gamma is a Chi- 
square type procedure which has inherent problems which 
does not lend itself to easy data independent analysis. 
Instead it requires judicious selection of parameters. 

This analysis did not take advantage of this option, hence, 
the rejection of the gamma could be attributable to poor 
cell location choices. 

Every data set was analyzed using a test for equality of 
shape parameters as proposed by Tubbs (1983) . This hypothe- 
sis of equal shape parameters was rejected in rvery case. 
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APPENDIX A 


Univariate summary statistics generated using 
PROC UNIVARIATE is the Statistical ^alysis 
System. 
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The next, four pages contain the sunsnary statistics 
for each of the univariate variables Au, Lu, Av, and Lv, 
respectively. The reference altitude is 

ALTITUDE - 4000 



4 
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The next four pages contain the summary statistics 
for each of the univariate variables Au, Lu, Av, and Lv, 


respectively. 
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The next four pages contain the summary 
for each of the univariate variables Au, Lu 
respectively. The reference altitude is 
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The next four pages contain the summary statistics 
for each of the univariate variables Au, Lu, Av, and Lv, 
respectively. The reference altitude is 
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The next four pages contain the sximinary statistics 
for each of the univariate variables Au, Lu, Av, and Lv, 
respectively. The reference altitude is 
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The next four pages contain the svamnary statistics 
for each of the xmivariate variables Au, Lu, Av, and Lv, 
respectively. The reference altitude is 
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APPENDIX B 


Plots for the bivariate gamma density function, scatter 
plots, and contours are given for each data set. 
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